Viewing STIS Data

0. Introduction

The Space Telescope Imaging Spectrograph (STIS) is a versatile imaging spectrograph installed on the Hubble Space Telescope (HST), covering a wide range of wavelengths from the near-infrared region into the ultraviolet.

This tutorial aims to prepare new users to begin analyzing STIS Data by going through downloading data, reading and viewing spectrum, and viewing STIS image.

There are three detectors on STIS: FUV-MAMA, NUV-MAMA, and CCD. While the detectors are designed for different scientific purposes and operate at different wavelength, their data are organized in the same structure. Thus we are using FUV-MAMA data as an example in this notebook, and leave the rest as exercises to the users.

For detailed overview of the STIS instrument, see the STIS Instrument Handbook. \ For more information on STIS data analysis and operations, see the STIS Data Handbook.

Defining some terms:

0.1 Import necessary packages

We will import the following packages:

  • astropy.io fits for accessing FITS files
  • astropy.table Table for creating tidy tables of the data
  • astroquery.mast Observations for finding and downloading data from the MAST archive
  • pathlib for managing system paths
  • matplotlib.pyplot for plotting data
  • IPython.display for formatting display
  • numpy to handle array functions
  • pandas to make basic tables and dataframes
  • stistools for quick operations on STIS Data
In [1]:
# Import for: Reading in fits file
from astropy.table import Table
from astropy.io import fits

# Import for: Downloading necessary files. (Not necessary if you choose to collect data from MAST)
from astroquery.mast import Observations

# Import for: Managing system variables and paths
from pathlib import Path

# Import for: Plotting and specifying plotting parameters
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.cm as cm
from IPython.display import display

# Import for: Quick Calculation and Data Analysis
import numpy as np
import pandas as pd
import math

# Import for operations on STIS Data
import stistools
The following tasks in the stistools package can be run with TEAL:
   basic2d      calstis     ocrreject     wavecal        x1d          x2d
/Users/kding/miniconda3/envs/stis/lib/python3.7/site-packages/stsci/tools/nmpfit.py:8: UserWarning: NMPFIT is deprecated - stsci.tools v 3.5 is the last version to contain it.
  warnings.warn("NMPFIT is deprecated - stsci.tools v 3.5 is the last version to contain it.")
/Users/kding/miniconda3/envs/stis/lib/python3.7/site-packages/stsci/tools/gfit.py:18: UserWarning: GFIT is deprecated - stsci.tools v 3.4.12 is the last version to contain it.Use astropy.modeling instead.
  warnings.warn("GFIT is deprecated - stsci.tools v 3.4.12 is the last version to contain it."

1. Downloading STIS data from MAST using astroquery

There are other ways to download data from MAST such as using CyberDuck. We are only showing how to use astroquery in this notebook

In [2]:
# make directory for downloading data
datadir = Path('./data')
datadir.mkdir(exist_ok=True)
In [3]:
# Search target objscy by obs_id
target = Observations.query_criteria(obs_id='ODJ102010')
# get a list of files assiciated with that target
FUV_list = Observations.get_product_list(target)
# Download only the SCIENCE fits files
Observations.download_products(FUV_list,productType="SCIENCE",extension='fits',download_dir=str(datadir))
INFO: Found cached file data/mastDownload/HST/odj102010/odj102010_raw.fits with expected size 8429760. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odj102010/odj102010_x1d.fits with expected size 77760. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odj102010/odj102010_x2d.fits with expected size 14474880. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odj102010/odj102010_flt.fits with expected size 10535040. [astroquery.query]
Out[3]:
Table length=4
Local PathStatusMessageURL
str50str8objectobject
data/mastDownload/HST/odj102010/odj102010_raw.fitsCOMPLETENoneNone
data/mastDownload/HST/odj102010/odj102010_x1d.fitsCOMPLETENoneNone
data/mastDownload/HST/odj102010/odj102010_x2d.fitsCOMPLETENoneNone
data/mastDownload/HST/odj102010/odj102010_flt.fitsCOMPLETENoneNone

2. Reading in the data

2.1 Investigating the data - Basics

Before doing any operation on the data, we want to first explore the basics and data file structures.

The Aperture extracted, background subtracted, flux and wavelength calibrated spectra data is stored in fits file with suffix _x1d. While we are using _x1d fits file as an example of investigating STIS table data, the following method of reading in data and viewing spectrum or other fields can be applied to any table data, either calibrated or uncalibrated. For more information on STIS file naming conventions, see Types of STIS Files.

Open the x1d fits file and explore its info and header:

In [4]:
#get information about the fits file
x1d_file=Path("./data/mastDownload/HST/odj102010/odj102010_x1d.fits")
fits.info(x1d_file)
Filename: data/mastDownload/HST/odj102010/odj102010_x1d.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     274   ()      
  1  SCI           1 BinTableHDU    156   1R x 19C   [1I, 1I, 1024D, 1024E, 1024E, 1024E, 1024E, 1024E, 1024E, 1024I, 1E, 1E, 1I, 1E, 1E, 1E, 1E, 1024E, 1E]   

The primary header that stores keyword information describing the global properties of all of the exposures in the file

In [5]:
#get header of the fits file
x1d_header_0 = fits.getheader(x1d_file,0)

try:
    print(x1d_header_0["instrument"])
except KeyError:
    print(x1d_header_0["primesi"])
for key in ["detector","obsmode", "opt_elem", "targname"]:
    print("{0}: {1}".format(key, x1d_header_0[key]))
print()   
STIS
detector: FUV-MAMA
obsmode: ACCUM
opt_elem: G140M
targname: -ALF-CEP

You can change the keys to check for other fields and metadata, or directly print the x1d_header to get all header information.

Some other metadata, such as exposure data and time, are stored in the first extension.

In [6]:
x1d_header_1 = fits.getheader(x1d_file,1)

date = x1d_header_1["DATE-OBS"]
time = x1d_header_1["Time-OBS"]
exptime = x1d_header_1["EXPTIME"]

print("The data were taken on {d}, starting at {t}, with the exposure time of {e} seconds".format(d=date,t=time,e=exptime))
The data were taken on 2018-01-14, starting at 08:31:01, with the exposure time of 900.0 seconds

2.2 Reading the table data

The main science data is stored at the first extension of the x1d fits file. We first read in the data as an astropy table.

In [7]:
#get data
x1d_data = Table.read(x1d_file,1)
# Display a representation of the data table:
x1d_data
WARNING: UnitsWarning: 'erg/s/cm**2/Angstrom' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Out[7]:
Table masked=True length=1
SPORDERNELEMWAVELENGTH [1024]GROSS [1024]BACKGROUND [1024]NET [1024]FLUX [1024]ERROR [1024]NET_ERROR [1024]DQ [1024]A2CENTEREXTRSIZEMAXSRCHBK1SIZEBK2SIZEBK1OFFSTBK2OFFSTEXTRLOCY [1024]OFFSET
AngstromsCounts/sCounts/sCounts/serg / (Angstrom cm2 s)erg / (Angstrom cm2 s)Counts/spixpixpixpixpixpixpixpixpix
int16int16float64float32float32float32float32float32float32int16float32float32int16float32float32float32float32float32float32
110241513.611692349868 .. 1567.6686700923991.423279 .. 0.048600590.0008550065 .. -6.92308e-051.422424 .. 0.048669821.96904e-12 .. 9.89624e-145.507516e-14 .. 1.495948e-140.03978601 .. 0.0073570892564 .. 2564389.577711102455-300300382.2857 .. 397.2364351.8867
In [8]:
# We can also get the columns of this table:
columns = x1d_data.colnames
columns
Out[8]:
['SPORDER',
 'NELEM',
 'WAVELENGTH',
 'GROSS',
 'BACKGROUND',
 'NET',
 'FLUX',
 'ERROR',
 'NET_ERROR',
 'DQ',
 'A2CENTER',
 'EXTRSIZE',
 'MAXSRCH',
 'BK1SIZE',
 'BK2SIZE',
 'BK1OFFST',
 'BK2OFFST',
 'EXTRLOCY',
 'OFFSET']

Another popular way of reading in fits data from hdu list as "FITS_rec":

In [9]:
with fits.open(x1d_file) as hdulist:
    fuv_x1d_data = hdulist[1].data

3. Plotting the spectrum

3.1 Making a simple plot of the spectrum

The actual data of each columns are stored in arrays with equal lengths. We collect the actual spectrum data and plot it, together with the error bar.

In [10]:
# From the astropy table, we first get all the data we need: wavelength, flux, and error
# notice that for astropy table, the column name is case sensitive
wl,flux,err=x1d_data[0]["WAVELENGTH","FLUX","ERROR"]

# make a plot of the data, use this cell to specify the format of the plot.
matplotlib.rcParams['figure.figsize'] = (20,7)
plt.style.use("bmh")

plt.plot(wl, flux, #the x-data, y-data, and y-axis error
             marker=".",markersize="2",markerfacecolor='w', markeredgewidth=0) #specifies the data points style
plt.title("STIS FUV Spectrum")
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")
Out[10]:
Text(0, 0.5, 'Flux [ergs/s/cm$^2$/Ã…]')

You can zoom in to a specific wavelength range using a mask

In [11]:
# create a mask on wavelength between 1540 Angstrom and 1560 Angstrom
zoom_mask = (wl > 1540) & (wl < 1560)
# apply the created mask to the wavelength, flux, and error
wl_region, flux_region, err_region = wl[zoom_mask],flux[zoom_mask],err[zoom_mask]
# plot the spectrum in that region
plt.plot(wl_region, flux_region, #the x-data, y-data, and y-axis error
             marker=".",markersize="2",markerfacecolor='w', markeredgewidth=0, #specifies the data points style
        color="blue") #specifies the format of lines
plt.title("STIS FUV Spectrum")
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")
Out[11]:
Text(0, 0.5, 'Flux [ergs/s/cm$^2$/Ã…]')

You can also plot the error bar together with the spectrum. For more errorbar styling options, see matplotlib.pyplot.errorbar

In [12]:
plt.errorbar(wl_region, flux_region, err_region,#the x-data, y-data, and y-axis error
             marker=".",markersize="2",markerfacecolor='w', markeredgewidth=0,color="blue", #specifies the data points style
             ecolor="dodgerblue",capsize=0.1) #specifies the format of lines and error bar
plt.title("STIS FUV Spectrum")
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")
Out[12]:
Text(0, 0.5, 'Flux [ergs/s/cm$^2$/Ã…]')

For more information on formatting the plots using matplotlib, see matplotlib.pyplot.plot, matplotlib.pyplot.errorbar.

4.Working with data quality flags

Data quality flags are assigned to each pixel in the data quality extension. Each flag has a true (set) or false (unset) state. Flagged conditions are set as specific bits in a 16-bit integer word. For a single pixel, this allows for up to 15 data quality conditions to be flagged simultaneously, using the bitwise logical OR operation. Note that the data quality flags cannot be interpreted simply as integers but must be converted to base-2 and interpreted as flags. These flags are set and used during the course of calibration, and may likewise be interpreted and used by downstream analysis applications.

4.1 Data quality frequencies histogram

Make a histogram according to the data quality flags, and label the bins by what each data quality values actually means. More info: https://hst-docs.stsci.edu/stisdhb/chapter-2-stis-data-structure/2-5-error-and-data-quality-array

In [13]:
#First get the data quality flag from the x1d fits file, and convert them to log2 values
#The data quality flag is a masked array that "hides" the pixels with no data quality issue. 
#We fill those "good points" with -1 in our case
x1d_dq = x1d_data[0]["DQ"].filled(0)
dq_bits = []
for dq in x1d_dq:
    if dq == 0:
        dq_bits.append(-1)
    else:
        for b in range(0,15):
            if dq & (2**b):
                dq_bits.append(b)
In [14]:
#Assign the meaning of each data quality and make the histogram
meanings = ["No Anomalies", "Error in the Reed Solomon decoding", "Lost data replaced by fill values",
           "Bad detector pixel (e.g., bad column or row, mixed science and bias for overscan, or beyond aperture)",
           "Data masked by occulting bar","Pixel having dark rate > 5 σ times the median dark level",
           "Large blemish, depth > 40% of the normalized p-flat (repeller wire)","Vignetted pixel",
           "Pixel in the overscan region",
           "Saturated pixel, count rate at 90% of max possible—local non-linearity turns over and is multi-valued; \
pixels within 10% of turnover and all pixels within 4 pixels of that pixel are flagged.",
           "Bad pixel in reference file","""Small blemish, depth between 40% and 70% of the normalized flat. \
Applies to MAMA and CCD p-flats""",
           ">30% of background pixels rejected by sigma-clip, or flagged, during 1-D spectral extraction",
           "Extracted flux affected by bad input data",
            "Data rejected in input pixel during image combination for cosmic ray rejection",
           "Extracted flux not CTI corrected because gross counts are ≤ 0"]

flags = [int(2**i) for i in range(-1,15)]
bins_ = np.arange(-1, 15)
plt.title("Data quality flag frequencies")
h = plt.hist(dq_bits,bins=bins_)
x = plt.xticks(bins_+0.5, labels=flags, rotation='vertical', fontsize=15)
plt.show()

bits = ['0'*(15-len(bin(i)[2:]))+bin(i)[2:] for i in flags]
dq_table = pd.DataFrame({"FLAG Value":flags,"Bit Setting":bits,"Quality Condition Indicated":meanings})

pd.set_option('display.max_colwidth', None)
display(dq_table)
FLAG Value Bit Setting Quality Condition Indicated
0 0 000000000000000 No Anomalies
1 1 000000000000001 Error in the Reed Solomon decoding
2 2 000000000000010 Lost data replaced by fill values
3 4 000000000000100 Bad detector pixel (e.g., bad column or row, mixed science and bias for overscan, or beyond aperture)
4 8 000000000001000 Data masked by occulting bar
5 16 000000000010000 Pixel having dark rate > 5 σ times the median dark level
6 32 000000000100000 Large blemish, depth > 40% of the normalized p-flat (repeller wire)
7 64 000000001000000 Vignetted pixel
8 128 000000010000000 Pixel in the overscan region
9 256 000000100000000 Saturated pixel, count rate at 90% of max possible—local non-linearity turns over and is multi-valued; pixels within 10% of turnover and all pixels within 4 pixels of that pixel are flagged.
10 512 000001000000000 Bad pixel in reference file
11 1024 000010000000000 Small blemish, depth between 40% and 70% of the normalized flat. Applies to MAMA and CCD p-flats
12 2048 000100000000000 >30% of background pixels rejected by sigma-clip, or flagged, during 1-D spectral extraction
13 4096 001000000000000 Extracted flux affected by bad input data
14 8192 010000000000000 Data rejected in input pixel during image combination for cosmic ray rejection
15 16384 100000000000000 Extracted flux not CTI corrected because gross counts are ≤ 0

4.2 Removing "Serious Data Quality"

Through the calibaration pipeline, some data qualities are marked "serious". The value of serious data qualities are marked through "SDQFLAGS". We can decompose that value according to the bits in order to see the specific data quality flags that are marked serious.

In [15]:
sdqFlags_fuv=fits.getheader(x1d_file,1)["SDQFLAGS"]
print("The SDQFLAGS is {flag}, which is in binary {bin_flag},".format(flag=sdqFlags_fuv,bin_flag=bin(sdqFlags_fuv)[2:]))
print("Therefore the following data qualities are marked \"serious\":")
for i in range(15):
    if 2**(i) & sdqFlags_fuv:
        print("\t"+str(i+1)+":" +meanings[i+1])
The SDQFLAGS is 31743, which is in binary 111101111111111,
Therefore the following data qualities are marked "serious":
	1:Error in the Reed Solomon decoding
	2:Lost data replaced by fill values
	3:Bad detector pixel (e.g., bad column or row, mixed science and bias for overscan, or beyond aperture)
	4:Data masked by occulting bar
	5:Pixel having dark rate > 5 σ times the median dark level
	6:Large blemish, depth > 40% of the normalized p-flat (repeller wire)
	7:Vignetted pixel
	8:Pixel in the overscan region
	9:Saturated pixel, count rate at 90% of max possible—local non-linearity turns over and is multi-valued; pixels within 10% of turnover and all pixels within 4 pixels of that pixel are flagged.
	10:Bad pixel in reference file
	12:>30% of background pixels rejected by sigma-clip, or flagged, during 1-D spectral extraction
	13:Extracted flux affected by bad input data
	14:Data rejected in input pixel during image combination for cosmic ray rejection
	15:Extracted flux not CTI corrected because gross counts are ≤ 0

We can now remove the data points with SDQ flags. In the plots below, the first one is the same as the spectrum in which we do not handle the SDQ flags. The second one is the spectrum with SDQ flagged data removed, and the SDQ flagged data points are marked with red '+'.

In [16]:
matplotlib.rcParams['figure.figsize'] = (20,10)
plt.style.use("bmh")

# First zoom in to the region where SDQ got removed
wl,flux,err,dq=x1d_data[0]["WAVELENGTH","FLUX","ERROR","DQ"]
zoom_mask = (wl > min(wl)) & (wl < 1520)
wl_zoom, flux_zoom, err_zoom, dq_zoom = wl[zoom_mask],flux[zoom_mask],err[zoom_mask],dq[zoom_mask]

# Filter the datapoints to where there are no serious DQ flags
mask_noSDQ = np.ones(len(dq_zoom),dtype=bool)
for i in range(0,len(dq_zoom)):
    if dq_zoom[i] & sdqFlags_fuv:
            mask_noSDQ[i]=False
# get the spectrum without SDQ using the mask we just created
wvln_noSDQ, flux_noSDQ,err_noSDQ = wl_zoom[mask_noSDQ], flux_zoom[mask_noSDQ],err_zoom[mask_noSDQ]
# inverse the _noSDQ mask to collect the data points with SDQ flags
mask_SDQ = [not elem for elem in mask_noSDQ]
wvln_SDQ, flux_SDQ = wl_zoom[mask_SDQ], flux_zoom[mask_SDQ]

# plot1: the spectrum with SDQ flagged data included
plt.subplot(2,1,1)
plt.plot(wl_zoom, flux_zoom, #the x-data, y-data, and y-axis error
             marker=".",markersize="2",markerfacecolor='w', markeredgewidth=0) #specifies the data points style
plt.title("FUV with SDQ flagged data")
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")

# plot2: plot the spectrum without SDQ flagged data, then mark the SDQ data points with +
plt.subplot(2,1,2)
# Plot the filtered datapoints
plt.plot(wvln_noSDQ, flux_noSDQ, #the x-data, y-data, and y-axis error
             marker=".",markersize="2",markerfacecolor='w', markeredgewidth=0,label="SDQ removed spectrum")
plt.plot(wvln_SDQ, flux_SDQ,'r+',label="SDQ flagged data")
plt.legend()

# Format the figure
plt.title("FUV without SDQ flagged data")
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")
plt.tight_layout()

5. Visualizing STIS Image

The STIS images are stored as two-dimensional arrays in FITS image extension files. For more information on STIS image files and extension, see STIS FITS Image Extension Files

5.1 Exploring image file structure

The rectified, wavelength and flux calibrated first order spectra or Geometrically corrected imaging data is stored in the fits file with the x2d extension. Similar to what we did to the x1d file, we first open the fits file to explore its file structure.

In [17]:
#read in the x2d file and get its info
x2d_file=Path("./data/mastDownload/HST/odj102010/odj102010_x2d.fits")
fits.info(x2d_file)
Filename: data/mastDownload/HST/odj102010/odj102010_x2d.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     277   ()      
  1  SCI           1 ImageHDU       120   (1201, 1201)   float32   
  2  ERR           1 ImageHDU        61   (1201, 1201)   float32   
  3  DQ            1 ImageHDU        44   (1201, 1201)   int16   
  • The first, of extension type SCI, stores the science values.
  • The second, of extension type ERR, contains the statistical errors, which are propagated through the calibration process. It is unpopulated in raw data files.
  • The third, of extension type DQ, stores the data quality values, which flag suspect pixels in the corresponding SCI data. It is unpopulated in raw data files.

ch2_stis_data4.1.jpg

Similarly, we can also get the header from this fits file to see the scientific metadata.

In [18]:
#get header of the fits file
x2d_header = fits.getheader(x2d_file,0)

5.2 Showing the image

Now we collect the science image data from the fits file and show the image.

In [19]:
# get data as a numpy array
with fits.open(x2d_file) as hdu_list:
    x2d_data = hdu_list[1].data

Make a histogram on the magnitude of the image data so that we have a general idea on the distribution. Knowing the distribution is essential for us to normalize the data when showing the image.

In [20]:
# remove infinities and nans from the data when making the histogram
filtered_data = [v for row in np.log10(x2d_data) for v in row if not (math.isinf(v) or math.isnan(v))]
plt.rcParams["figure.figsize"] = plt.rcParamsDefault["figure.figsize"]
plt.hist(filtered_data)
plt.xlabel("Magnitude order of _x2d image data")
plt.ylabel("count")
plt.title("_x2d Data Magnitude Order frequencies")
/Users/kding/miniconda3/envs/stis/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in log10
  
/Users/kding/miniconda3/envs/stis/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in log10
  
Out[20]:
Text(0.5, 1.0, '_x2d Data Magnitude Order frequencies')

When showing the image, we normalize the color of each pixel to a specific range through vmin and vmax to make the features of image clear. These values typically matches the magnitude of the x2d data according to the histogram above, but can be experimented and changed to bring out features at different levels.

In [21]:
# show the image
matplotlib.rcParams['figure.figsize'] = (10,10)
plt.imshow(x2d_data,origin='lower',vmin=0,vmax=1e-13,cmap="viridis")
Out[21]:
<matplotlib.image.AxesImage at 0x7f986a0aedd0>

For more color map and normalizing options, see: Choosing Colormaps in Matplotlib, matplotlib.pyplot.imshow.

5.3 Removing serious data quality pixels

In [22]:
# get the serious data quality flag
sdqFlags_fuv=fits.getheader(x2d_file,1)["SDQFLAGS"]
# get data quality flags of each pixels
with fits.open(x2d_file) as hdu_list:
    x2d_dq = hdu_list[3].data
# create a mask of bad pixels and set them to nan
def check_dq(dq):
    return bool(dq & sdqFlags_fuv)
mask = np.vectorize(check_dq)(x2d_dq)
x2d_mask = np.ma.array(x2d_data,mask=mask,fill_value=np.nan)
# plot the image
plt.imshow(x2d_mask,origin='lower',vmin=0,vmax=1e-13,cmap="viridis")
Out[22]:
<matplotlib.image.AxesImage at 0x7f9888c5f250>

6.Working with Time-Tag data

The MAMA detecters have a unique Time-Tag mode besides ACCUM mode. TIME-TAG mode is used for high-time-resolution spectroscopy and imaging in the UV. In TIME-TAG mode, the position and detection time of every photon is recorded in an event list. The Time-Tag mode operation for the MAMA detectors can be found at: MAMA TIME-TAG Mode.

In TIME-TAG mode, the position and detection time of every photon is recorded in an event list. First collect the _tag data:

In [23]:
# Search target objscy by obs_id
target = Observations.query_criteria(obs_id='odgxt9010')
# get a list of files assiciated with that target
NUV_list = Observations.get_product_list(target)
# Download only the SCIENCE fits files
Observations.download_products(NUV_list,extension='fits',download_dir=str(datadir))
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_jif.fits with expected size 34560. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_jit.fits with expected size 103680. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_jwf.fits with expected size 25920. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_jwt.fits with expected size 20160. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_spt.fits with expected size 241920. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_tag.fits with expected size 48064320. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_trl.fits with expected size 37440. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_wav.fits with expected size 4262400. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_wsp.fits with expected size 112320. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_asn.fits with expected size 11520. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_raw.fits with expected size 16819200. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_x1d.fits with expected size 77760. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_x2d.fits with expected size 14474880. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/odgxt9010/odgxt9010_flt.fits with expected size 10535040. [astroquery.query]
Out[23]:
Table length=14
Local PathStatusMessageURL
str50str8objectobject
data/mastDownload/HST/odgxt9010/odgxt9010_jif.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_jit.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_jwf.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_jwt.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_spt.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_tag.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_trl.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_wav.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_wsp.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_asn.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_raw.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_x1d.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_x2d.fitsCOMPLETENoneNone
data/mastDownload/HST/odgxt9010/odgxt9010_flt.fitsCOMPLETENoneNone

6.1 Investigating the _tag data

In [24]:
# get info about the tag extension fits file
tag = Path("./data/mastDownload/HST/odgxt9010/odgxt9010_tag.fits")
fits.info(tag)
Filename: data/mastDownload/HST/odgxt9010/odgxt9010_tag.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     218   ()      
  1  EVENTS        1 BinTableHDU    147   4802131R x 4C   [1J, 1I, 1I, 1I]   
  2  GTI           1 BinTableHDU     22   1R x 2C   [1D, 1D]   

The _tag fits file has two binary table extensions: EVENTS and GTI.

In [25]:
# get header of the EVENTS extension
# print only the TIMETAG EVENTS TABLE COLUMNS (line 130-147)
fits.getheader(tag,1)[130:147]
Out[25]:
                                                                                
              / TIMETAG EVENTS TABLE COLUMNS                                    
                                                                                
TTYPE1  = 'TIME    '           / event clock time                               
TFORM1  = '1J      '           / data format for TIME: 32-bit integer           
TUNIT1  = 'seconds '           / units for TIME: seconds                        
TSCAL1  =             0.000125 / scale factor to convert s/c clock ticks to sec 
TZERO1  =                  0.0 / TIME zero point: starting time of obs.         
TTYPE2  = 'AXIS1   '           / Doppler corrected axis 1 detector coordinate   
TFORM2  = '1I      '           / data format for AXIS1: 16-bit integer          
TUNIT2  = 'pixels  '           / physical units for AXIS1: pixels               
TTYPE3  = 'AXIS2   '           / axis 2 detector coordinate                     
TFORM3  = '1I      '           / data format for AXIS2: 16-bit integer          
TUNIT3  = 'pixels  '           / physical units for AXIS2: pixels               
TTYPE4  = 'DETAXIS1'           / raw 1 detector coordinate                      
TFORM4  = '1I      '           / data format for DETAXIS1: 16-bit integer       
TUNIT4  = 'pixels  '           / physical units for DETAXIS1: pixels            

Columns in the EVENTS extension:

  • TIME: the time each event was recorded relative to the start time
  • AXIS1: pixel coordinate along the spectral axis with the corretion term on Doppler shifts
  • AXIS2: pixel coordinate along the spatial axis
  • DETAXIS1: pixel coordinate along the spectral axis without the corretion term on Doppler shifts
In [26]:
# get header of the GTI extension
fits.getheader(tag,2)
Out[26]:
XTENSION= 'BINTABLE'           / extension type                                 
BITPIX  =                    8 / bits per data value                            
NAXIS   =                    2 / number of data axes                            
NAXIS1  =                   16 / length of first data axis                      
NAXIS2  =                    1 / length of second data axis                     
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    2 / number of fields in each table row             
INHERIT =                    T / inherit the primary header                     
EXTNAME = 'GTI     '           / extension name                                 
EXTVER  =                    1 / extension version number                       
ROOTNAME= 'odgxt9010                         ' / rootname of the observation set
EXPNAME = 'odgxt9rfq                ' / exposure identifier                     
                                                                                
              / TIMETAG GTI TABLE COLUMN DESCRIPTORS                            
                                                                                
TTYPE1  = 'START   '           / start of good time interval                    
TFORM1  = '1D      '           / data format for START: REAL*8                  
TUNIT1  = 'seconds '           / physical units for START                       
TTYPE2  = 'STOP    '           / end of good time interval                      
TFORM2  = '1D      '           / data format for STOP: REAL*8                   
TUNIT2  = 'seconds '           / physical units for STOP                        
In [27]:
Table.read(tag,2)
Out[27]:
Table length=1
STARTSTOP
secondsseconds
float64float64
0.0226252380.222

Columns in the GTI extension:

  • START: start of good time interval
  • STOP: end of good time interval

Now we make a time series plot of the total flux over all wavelengths to see how the flux changes over the time interval:

In [28]:
# read the events data as a pandas dataframe for easier manipulation
events = Table.read(tag,1).to_pandas()
# get the good time interval from the GTI extension
start,stop = Table.read(tag,2)["START"],Table.read(tag,2)["STOP"]
# group the events by time with bin = 3 seconds
time = np.arange(int(start),int(stop),3)
ind = np.digitize(events['TIME'], time)
total_flux = events.groupby(ind).count()["TIME"]
# plot the flux as a function of time
# notice here that the unit of flux is counts since we are counting the number of events in a time series
matplotlib.rcParams['figure.figsize'] = (20,7)
plt.plot(time,total_flux,marker=".",mfc="w")
plt.xlabel("Time [s]")
plt.ylabel("Total Flux [counts]")
Out[28]:
Text(0, 0.5, 'Total Flux [counts]')

As the plot shows, though the total flux fluctuates, it is roughly a constant over the good time interval.

6.2 Converting Time_Tag into ACCUM image

Time tag data can be converted into ACCUM image using the inttag method in stistools. More information: Error and Data Quality Array

In [29]:
# define the output file directory
accum = "./data/mastDownload/HST/odgxt9010/odgxt9010_accum.fits"
# convert Time_Tag into ACCUM
# the first parameter is the path to the _tag fits file, the second parameter is the output directory
stistools.inttag.inttag(tag, accum)
imset: 1, start: 0.022625, stop: 2380.222, exposure time: 2380.199375

Then the output file is in the same structure as a STIS image fits file, which we can open and plot in the same way we explored above:

In [30]:
with fits.open(accum) as hdul:
    im = hdul[1].data
plt.imshow(im,origin='lower',vmin=0,vmax=6,cmap="viridis")
Out[30]:
<matplotlib.image.AxesImage at 0x7f98611aff90>

inttag with multiple data sets: rcount specifies the number of imsets, imcrements specifies the time interval for each imsets in seconds

In [31]:
stistools.inttag.inttag(tag,accum, rcount = 3, increment = 700)
fits.info(accum)
imset: 1, start: 0.022625, stop: 700.022625, exposure time: 700.0
imset: 2, start: 700.022625, stop: 1400.022625, exposure time: 700.0000000000001
imset: 3, start: 1400.022625, stop: 2100.022625, exposure time: 700.0
Filename: ./data/mastDownload/HST/odgxt9010/odgxt9010_accum.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     218   ()      
  1  SCI           1 ImageHDU       138   (1024, 1024)   float64   
  2  ERR           1 ImageHDU       138   (1024, 1024)   float64   
  3  DQ            1 ImageHDU       139   ()      
  4  SCI           2 ImageHDU       138   (1024, 1024)   float64   
  5  ERR           2 ImageHDU       138   (1024, 1024)   float64   
  6  DQ            2 ImageHDU       139   ()      
  7  SCI           3 ImageHDU       138   (1024, 1024)   float64   
  8  ERR           3 ImageHDU       138   (1024, 1024)   float64   
  9  DQ            3 ImageHDU       139   ()      

Compare the 3 accum images produced by inttag, using the same scale and colormap:

In [32]:
for i in range(1,8,3):
    with fits.open(accum) as hdul:
        im = hdul[i].data
    plt.subplot(1,3,int(i/3)+1)
    plt.imshow(im,origin='lower',vmin=0,vmax=2,cmap="viridis")
    plt.title("extension {}".format(i))
plt.tight_layout()

The output file is a series of extensions with each imset having a SCI, ERR, and DQ extension, as shown above.

7.Working with STIS Echelles data

The STIS Echelles mode uses diffraction and interference to separate a spectrum into different spectral orders (the keyword in the headers is SPORDER), with each spectral order covering different wavelength range. The echelles were designed to maximize the spectral range covered in a single echellogram. For more information, see Long-Slit Echelle Spectroscopy.

We first collect an Echelles data using astroquery:

In [33]:
# Search target objscy by obs_id
target = Observations.query_criteria(obs_id='OCTX01030')
# get a list of files assiciated with that target
NUV_list = Observations.get_product_list(target)
# Download only the SCIENCE fits files
Observations.download_products(NUV_list,extension='fits',download_dir=str(datadir))
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_jif.fits with expected size 216000. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_jit.fits with expected size 164160. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_jwf.fits with expected size 17280. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_jwt.fits with expected size 11520. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_sfl.fits with expected size 10535040. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_spt.fits with expected size 371520. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_trl.fits with expected size 43200. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_wav.fits with expected size 2139840. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_wsp.fits with expected size 69120. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_asn.fits with expected size 11520. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_raw.fits with expected size 67317120. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_x1d.fits with expected size 7329600. [astroquery.query]
INFO: Found cached file data/mastDownload/HST/octx01030/octx01030_flt.fits with expected size 84139200. [astroquery.query]
Out[33]:
Table length=13
Local PathStatusMessageURL
str50str8objectobject
data/mastDownload/HST/octx01030/octx01030_jif.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_jit.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_jwf.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_jwt.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_sfl.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_spt.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_trl.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_wav.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_wsp.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_asn.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_raw.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_x1d.fitsCOMPLETENoneNone
data/mastDownload/HST/octx01030/octx01030_flt.fitsCOMPLETENoneNone

7.1 Showing Echelles Image

Open the _flt (Flat-fielded science) image in the same way we open other image files, and show the image:

In [34]:
flt_file="./data/mastDownload/HST/octx01030/octx01030_flt.fits"
with fits.open(flt_file) as hdu:
        image = hdu[1].data
        plt.imshow(image,origin='lower',vmin=0,vmax=1,cmap="viridis")

As shown in the image above, there are multiple horizontal bands with different spectral orders. Each spectral order also has distinct wavelength range, which we will explore in the next section.

7.2 Plotting Echelles Spectrum

We first read in the _x1d data as an astropy table. Notice that when we read in the FUV _x1d data in section2.2, the table has a single row with SPORDER = 1. But for echelles mode data, the data is separated into multiple rows, with each row having a distinct order.

In [35]:
echelles_file="./data/mastDownload/HST/octx01030/octx01030_x1d.fits"
echelles_data = Table.read(echelles_file,1)
echelles_data
WARNING: UnitsWarning: 'erg/s/cm**2/Angstrom' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
Out[35]:
Table masked=True length=23
SPORDERNELEMWAVELENGTH [1024]GROSS [1024]BACKGROUND [1024]NET [1024]FLUX [1024]ERROR [1024]NET_ERROR [1024]DQ [1024]A2CENTEREXTRSIZEMAXSRCHBK1SIZEBK2SIZEBK1OFFSTBK2OFFSTEXTRLOCY [1024]OFFSET
AngstromsCounts/sCounts/sCounts/serg / (Angstrom cm2 s)erg / (Angstrom cm2 s)Counts/spixpixpixpixpixpixpixpixpix
int16int16float64float32float32float32float32float32float32int16float32float32int16float32float32float32float32float32float32
6710243021.857427057288 .. 3072.451049005796-0.0002917696 .. 0.83085710.003859486 .. 0.001177549-0.004151256 .. 0.8296795-8.482215e-15 .. 3.38419e-122.199675e-17 .. 4.074204e-131.076537e-05 .. 0.099884572628 .. 258071.193337105527.28738-28.0199575.95515 .. 66.737220.6031896
6810242977.380919809147 .. 3027.257338502545-0.008585975 .. 0.41904820.01126011 .. 0.002705485-0.01984609 .. 0.4163427-3.480781e-14 .. 1.039089e-121.027387e-16 .. 1.776893e-135.857771e-05 .. 0.071196642580 .. 2564126.36837105526.57536-27.28738131.2862 .. 122.41060.6031896
6910242934.193937243135 .. 2983.372006818285-0.0003262974 .. 0.97402620.01237513 .. 0.002681673-0.01270143 .. 0.9713445-2.004744e-14 .. 2.056522e-121.795185e-17 .. 2.19945e-131.137373e-05 .. 0.10388532564 .. 2564179.90447105525.88353-26.57536184.4561 .. 176.38990.5168002
7010242892.241202288543 .. 2940.739046013224-0.0003810195 .. 1.4136890.009411922 .. 0.004319549-0.009792942 .. 1.40937-1.405423e-14 .. 2.825297e-121.764158e-17 .. 2.479927e-131.229259e-05 .. 0.12370862564 .. 2564232.0847105525.21152-25.88353236.3518 .. 228.98390.5693123
7110242851.470552660116 .. 2899.305600220004-0.0002589348 .. 1.3539630.002474391 .. 0.005335689-0.002733326 .. 1.348627-3.772761e-15 .. 2.665274e-121.400337e-17 .. 2.393677e-131.014529e-05 .. 0.12111992564 .. 2564282.73527105524.55895-25.21152286.8198 .. 280.10420.5335968
7210242811.832724517215 .. 2859.021746954763-0.0002629836 .. 0.37120380.004915517 .. 0.001902789-0.0051785 .. 0.369301-7.105272e-15 .. 7.38581e-131.412153e-17 .. 1.340935e-131.029212e-05 .. 0.067048652564 .. 2564332.08117105523.92547-24.55895335.8839 .. 329.92760.579392
.........................................................
8410242409.809495426899 .. 2450.35395445196-0.009019177 .. -0.024645230.0007367581 .. 0.0002889372-0.009755936 .. -0.02493417-2.288945e-14 .. -1.151903e-131.405606e-16 .. 8.281269e-145.990969e-05 .. 0.017925692580 .. 2564832.81127105517.6788-18.11187833.6877 .. 835.05370.6772078
8510242381.432648611571 .. 2421.501716113727-0.0003795348 .. 0.012159110.000633091 .. 0.0002293428-0.001012626 .. 0.01192976-2.612651e-15 .. 6.205271e-143.15962e-17 .. 1.380587e-131.224623e-05 .. 0.026542082564 .. 2564868.10687105517.26004-17.6788868.6304 .. 870.57980.5922183
8610242353.715896795153 .. 2393.319815723954-0.001006938 .. 0.00031683020.0005307035 .. 0.0001866426-0.001537641 .. 0.0001301876-4.130699e-15 .. 7.576767e-165.358452e-17 .. 1.41783e-131.994669e-05 .. 0.024361822564 .. 2564902.69197105516.85524-17.26004903.1175 .. 905.54280.6101555
8710242326.636473828409 .. 2365.785161512079-0.0004271266 .. -0.036412350.0004274778 .. 0.000170473-0.0008546044 .. -0.03658282-2.335252e-15 .. -2.532413e-133.526203e-17 .. 1.021887e-131.290443e-05 .. 0.014762012628 .. 2564936.48247105516.46401-16.85524936.36 .. 939.48320.6185523
8810242300.172648559522 .. 2338.875710356606-0.0003803358 .. -0.026999660.0003694183 .. 0.000119634-0.0007497541 .. -0.0271193-1.964272e-15 .. -2.173019e-133.219402e-17 .. 1.381312e-131.228832e-05 .. 0.017238782628 .. 2564969.4647105516.08601-16.46401969.1404 .. 972.77430.5765846
8910242274.303666679877 .. 2312.570408922325-0.0004408245 .. -0.035755230.0001131443 .. 8.110702e-05-0.0005539688 .. -0.03583634-1.748883e-15 .. -3.811664e-134.183209e-17 .. 1.557065e-131.325056e-05 .. 0.014639142628 .. 25641001.7797105515.72085-16.086011001.03 .. 1005.4380.601592

Now we can plot the spectrum of all spectral orders in one plot, with each spectral order having a distinct color:

In [36]:
# plot the spectrum; the color of each SPORDER is specified through a matplotlib built-in colormap
matplotlib.rcParams['figure.figsize'] = (20,7)

for i in range(len(echelles_data)):
    plt.plot(echelles_data[i]['WAVELENGTH'],echelles_data[i]['FLUX'],color=cm.get_cmap('prism')(i/len(echelles_data)),label=echelles_data[i]['SPORDER'],alpha=0.7)
plt.xlabel('Wavelength [' + chr(197) +']')
plt.ylabel("Flux [ergs/s/cm$^2$/" + chr(197) +"]")
plt.legend(loc='best')
plt.title("STIS Echelles Mode Spectrum")
Out[36]:
Text(0.5, 1.0, 'STIS Echelles Mode Spectrum')

As the spectrum illustrates, each spectral order covers a specific wavelength range. Notice that the wavelength ranges could overlap with each other.


About this Notebook

Author: Keyi Ding

Updated On: 2022-07-18

This tutorial was generated to be in compliance with the STScI style guides and would like to cite the Jupyter guide in particular.

Citations

If you use astropy, matplotlib, astroquery, or numpy for published research, please cite the authors. Follow these links for more information about citations:


Top of Page Space Telescope Logo

In [ ]: